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Abstract 

We present molecular dynamics results for a two component, two-dimensional Lennard- Jones 
supercooled liquid near the glass transition. We find that the supercooled liquid is spatially 
heterogeneous and that there are long-lived clusters whose size distribution satisfies a scaling 
relation up to a cutoff. The similarity of several properties of the supercooled liquid to those 
of a mean-field glass-forming fluid near the spinodal suggests that the glass transition in the 
supercooled liquid is associated with an underlying thermodynamic instability. 

1. Introduction 

The characterization of supercooled liquids near the glass transition is an active area of 
research JT|. Outstanding unsolved problems include the possible existence of an underly- 
ing thermodynamic glass transition, the history dependence of the glass properties, and the 
mechanisms responsible for the large increase in relaxation times as the glass transition is 
approached. In this article we discuss our molecular dynamics simulations of a two compo- 
nent, two-dimensional (d = 2) Lennard- Jones supercooled liquid. A summary of some of our 
results has been published earlier Qj. 

To help the reader understand our interpretation of the molecular dynamics data and 
the questions we pose, we first review the behavior of a mean-field model of a structural 
glass transition ||. In the mean- field model, particles interact via a repulsive, two-body 
potential of the form Q V(r) = / y d <j)('-fr), where r is the interparticle distance and d is the 
spatial dimension. The range of the interaction is R = 7 _1 . In the limit R — > oo, the 
static properties of the uniform fluid are described exactly by mean- field theory II. The 
equilibrium structure function S(k) can be shown to satisfy the Ornstein-Zernicke form, 
S(k) = 1/[1 + Ppcp(k)}, where p is the particle density, = k B T, and 4>{k) is the Fourier 
transform of 4>(x). We choose the step potential, <f)(x) = 1 for x < 1 and <f>(x) = for x > 1. 
The Fourier transform (j){k) is negative for various values of k. This property of 4>{k) implies 
that in the mean-field limit, S(k) becomes negative in some range of k for sufficiently large 
(3 p. Hence for fixed p, the system has a spinodal singularity § at a temperature T = T s 
which is defined by the condition, 1 + /3 s p<ft(ko) = 0, where ko is the location of the minimum 
of 4>{k). At T = T s , S(k) diverges at k ; for T > T s , S(k ) diverges as (T — T s )~ 7 with the 
mean-field exponent 7 = 1. Note that the spinodal is the limit of thermodynamic stability 
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of the uniform phase and is a critical point (for a given density). However, unlike the usual 
critical point, the structure function diverges at a nonzero value of k. 

For d = 3 and p = 1.95, the spinodal is at T s = 0.705 {k B = 1) ||. Although the 
singularity is well defined only in the mean- field limit, simulations must be done for finite 
R. Monte Carlo simulations for R = 3 with p = 1.95 yield a S(k) which has a maximum at 
k 7^ that increases rapidly as T is decreased until T « 0.75, below which the peak ceases 
to increase as T is lowered [B[H|. This behavior is characteristic of a pseudospinodal [M. 

One way of understanding the simulation results for S(k) is to interpret the spinodal as 
a singularity of the free energy as a function of T, p. If the interaction range is finite, the 
singularity is in the complex (T, p) plane as has been found in transfer matrix studies of 
long-range ferromagnetic Ising models |§. The singularity moves closer to the real axis as 
the range of interaction increases. We refer to the complex singularity in finite range systems 
as a pseudospinodal. The physical effects of the pseudospinodal depend on how far it is from 
the real axis. As R increases, the pseudospinodal better approximates a true singularity and 
has measurable effects if R is sufficiently large. 

The fact that the mean-field model has a spinodal is unambiguous. We now discuss the 
reasons why we can associate the spinodal with a thermodynamic glass transition. In the 
Monte Carlo simulation for fixed p we equilibrate the system at a high T > T s where the 
uniform phase is stable and then quench it to T < T s where the uniform phase is unstable. 
After the quench the particles immediately form clumps with order pR? particles in each 
clump ||. That is, the system breaks up into regions of high particle density surrounded 
by regions of low particle density. The arrangement of the clumps is noncrystalline and 
their number depends on the quench history pM. The free energy has been calculated 
numerically in the mean-field limit and has many minima corresponding to different numbers 
of clumps [0,0. For T < T s , the system remains trapped in a local free energy minimum. 
The minimum the system chooses appears to depend on the quench history. 

The multi-minima free energy and the quench history dependence suggest that the mean- 
field model has a metastable glass phase for T <T S . However, the properties of a glass are 
associated more with its dynamical properties. As we will discuss, the glassy dynamics of 
the mean-field model is associated with the behavior of the clumps; the particles are not 
localized in the metastable glass phase. That is, the dynamical properties associated with 
single particle behavior do not show the usual signature of the approach to a glass. A simple 
argument based on the fact that all potential barriers in the mean-field model are finite || 
implies that the self-diffusion coefficient D is nonzero for all T > 0. Hence, if the observation 
time is sufficiently long, the mean square displacement of the particles increases indefinitely. 
The argument proceeds as follows. Because the separation between the clumps is order R, 
the interaction of the particles in a particular clump with the particles in neighboring clumps 
is minimized. Inside a clump, a particle undergoes a restricted random walk. A particle that 
attempts to leave a clump experiences a potential barrier due to the proximity of other 
clumps. Such a particle must interact with particles in the same clump and with particles 
in at least one other clump at some time during its possible escape. The upper bound of 
the potential barrier is cj R , where the constant c depends on d. (Recall that 7 d is the 
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strength of the interaction and R d is proportional to the number of particles in a clump such 
that 7^/?^ = 1.) The probability of leaving a clump in a Monte Carlo simulation is bounded 
from below by |e _C//fcsT for all R and hence D > 0. Similar arguments hold for a molecular 
dynamics simulation of the same system. 

The same argument implies that the mean-field model is ergodic for all T > if sin- 
gle particle properties are probed. The ergodic behavior can be characterized by several 
fluctuation metrics [ID]. The single particle energy fluctuation metric Q € (t) is given by 
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where 



e i (t) = U t e i (t')dt' (2) 
t Jo 

is the energy of particle i averaged over the time interval t, and (e(t)) = (1/N) E*Li ^(*)- 
The single particle energy q of particle i includes its kinetic energy and one-half of the 
potential energy of its interaction with other particles in the system. If the system is ergodic, 
Q e (t) ~ \jt for t sufficiently large lQfl . We find that Q e (t) exhibits ergodic behavior at 



T = 0.4, a value of T < T s . However, for T < 0.15, l/fi e (t) is not proportional to t during 
our longest runs, and the measured D is indistinguishable from zero at T = 0.15. Given 
our theoretical argument that D is nonzero for all T > 0, the observed behavior of l/Q e (t) 
for T < 0.15 implies only that the time for a particle to leave a clump is much longer than 
our observation time. We interpret this change from ergodic to nonergodic behavior as an 
apparent kinetic transition. 

How can we reconcile the T-dependence of D and Q e {t) with our identification of T s as 
the spinodal/glass transition? The answer lies in the dynamical properties of the clumps. 
For example, the diffusion coefficient of the center of mass of the clumps, D c \, is zero for 
T < T s in the mean-field limit. To understand this behavior, note that n c i, the mean number 
of particles in a clump, diverges as R d as R — > oo. In this limit, the number of clumps does 
not change with time and the mean numbers of particles exiting and entering a clump are 
equal. From the central limit theorem, the relative fluctuations of these quantities goes to 
zero as R and n d — > oo. We conclude that D c \ < D/^/n^., and hence D c \ = in the mean- 
field limit. Our simulations of D c \ for finite R are consistent with this prediction, and also 
indicate that the clumps do not see the same local environment, that is, they have different 
numbers of nearest neighbor clumps. This difference in the local environment persists as 
R —>■ oo because the clumps do not diffuse and cannot sample different local environments. 
Hence, the system is nonergodic on a clump (mass ~ R d ) scale for T < T s in the mean-field 
limit. 

In summary, the static properties of the mean-field glass phase are dominated by local- 
ized, long-lived structures (clumps) for T < T s , even though particles move from clump to 
clump. The time scale for the motion of the clumps diverges in the mean-field limit, and the 
spinodal/glass transition is seen dynamically only on a clump scale. However, we observe an 
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apparent kinetic transition that is associated with the slow diffusion of the particles and the 
finite duration of our runs. The temperature of this apparent transition is less than T s and 
depends on the observation time. 

The well characterized behavior of the mean-field glass model motivates us to ask if 
similar behavior occurs in systems with more realistic interactions. In Section 2 we discuss 
the static properties of a two component, two-dimensional system of Lennard- Jones particles 
and find that the system exhibits the usual properties associated with other simulations of 
glassy systems. In Section 3 we compute the static structure function S(k) and find evidence 
for a weak divergence in the height of the diffraction peak as T is lowered. This divergence 
is interpreted as evidence for the influence of a pseudospinodal. In Section 4 we propose 
a criterion for clusters of solid-like particles and find evidence of cluster scaling and an 
ergodic to nonergodic transition in the dynamics of the clusters. In Section 5 we discuss the 
interpretation of our molecular dynamics results in terms of a mean-field perspective. 

2. Lennard-Jones model. 

Because it is easier to identify geometrical structures in two dimensions than in three, we 
consider a two-dimensional system. We also specify that the particles interact via a Lennard- 
Jones potential because of the extensive simulations of supercooled Lennard-Jones liquids. 
However, a single component, two-dimensional Lennard-Jones system quickly nucleates to 



a solid and is unlikely to form a glass-like state. We follow Wong and Chester [ 1 1 1 who 
have done a Monte Carlo study of the quenched states of two component, two-dimensional 
Lennard-Jones systems with the goal of choosing the density and the size of the minority 
component so that it inhibits nucleation of the majority component to a solid. 

We designate by V ab the interaction between a particle belonging to species a and a 
particle belonging to species b and write V ab = Ae[(a ab /r) 12 — (a ab /r) 6 ]. We follow ref. |l 1 



and choose a aa = a, a bb = 3a/2, and a ab = \(a aa + a bb ) = 5a/ 4. The energy parameter 
e and the mass m are the same for both species. The Lennard-Jones potential is cutoff at 
r = 2.5 a and shifted so that the potential goes to zero at the cutoff. We choose units such 
that lengths are measured in terms of a, energies in terms of e, and the time in terms of 
r = (ma 2 /e) 1 / 2 . For example, the reduced number density is p* = pa 2 . In the following, 
we will omit the asterisk and quote all results in reduced units. Because the time can be 
expressed either in terms of the number of time steps or in terms of r, we will explicitly 
specify the units of time. 

The majority a component is taken to be 80% of the total number of particles N. All 
but one of the simulations are for N = 500. This relatively small number of particles was 
chosen so that long runs could be made. The duration of the run for each temperature was 
100 000 r unless otherwise noted. The molecular dynamics simulations were performed with 



periodic boundary conditions at constant energy and volume [12] at a density of p = 0.952 



corresponding to L = 22.92, where L is the linear dimension of the simulation cell. We 



used the velocity form of the Verlet algorithm [V2| with a time step of At = 0.005 r for all 
temperatures except for T = 5.5 for which we choose At = 0.0025 r to achieve reasonable 
energy conservation. The mean temperatures and pressures of our eight different runs are 
shown in Table |. As a typical example of the statistical accuracy of these mean values, 
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we note that the run at T = 3.1 corresponds to a mean temperature of 3.15 over the first 
50 000 r and 3.13 over the second 50 000 r. However for T = 2.25, the mean temperatures 
were 2.21 and 2.29 respectively, and it is possible that the lowest temperature run has not 
reached thermal equilibrium. 

Because our results are for constant density rather than constant pressure, the values of 
various quantities to be discussed here will differ slightly from our earlier results cited in 
ref. [Q] which were for a pressure of P ~ 70 and different densities. The initial configuration 
of each run at a particular temperature was the final configuration obtained previously [[ 



In some cases the positions had to be rescaled to obtain the desired density. Unless other- 
wise stated, the system was allowed to equilibrate for at least 50 000 r and then data was 
collected for 100 000 r. This duration is one to two orders of magnitude longer than reported 
previously |2[ and is equivalent to runs of approximately 2 x 10 -7 s for liquid Argon. 

We first compute the radial distribution function g(r) and the self-diffusion coefficient D 
and show that they exhibit behavior similar to that found in other simulations of supercooled 
liquids. All of the following results are for the majority particles only unless otherwise noted. 
The positions of the particles were saved every 5r and g{r) was computed in a separate 
program. In Fig. [l| we show g{r) for T = 2.7. Note the split second peak, a possible 
signature of a glassy state |TTJ; the split second peak of g{r) is not observed for T > 3.1. 
The increase in the height of the first peak of g(r) as the temperature is lowered (see Table |) 
indicates the growth of short-range order. We discuss the temperature dependence of the 
enthalpy and the constant pressure heat capacity in Appendix A. 

The single particle self-diffusion coefficient D is related to the mean square displacement 
R 2 {t) by lim f i? 2 (t) = 2dDt. An alternative way of determining D is from the velocity 



fluctuation metric VL v (t) defined as |T0 



N i d 
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= ^E^EEM*) - (Ut))} 2 , (3) 

i=l a=l 

where Uj jQ is the a-component of the velocity of particle i. The time average, Vi >a (t), of 
v it a is defined as in (fj) and the quantity (v a (t)) is the average of v it0l (t) over all particles. In 
our simulations the total momentum of the system is equal to zero, and hence (v a (t)) = 0. 
Because /„* Vj(t') df = r^t) - n(0), we can express Q v (t) as 

Q v (t) = R 2 {t)/2t 2 . (4) 

If the particles are undergoing diffusion, we have 

2D 

-> — • (5) 

The relation (||) was also used to determine D. 

For each value of T, the quantity |rj(t + t ) — rj(t )| 2 was computed in a separate program 
by grouping the particle coordinate data into blocks of duration 500 r and averaging over all 
possible choices of the origin to and combinations of time differences t < 250 for a particular 
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block. The results were then averaged over the 200 blocks. The velocity metric was computed 
at each time step in the main molecular dynamics program and the results also were averaged 
over 200 time origins. The typical time-dependence of Q v (0)/Q v (t) is shown in Fig. |2| for 
T = 2.25. Note that at this temperature 1/Q v (t) is not linear for t < 200 r, but becomes 
linear for longer observation times. 

Our results for D obtained from R 2 (t) and 1/Q v (t) are summarized in Table |I|. The 
difference in the values of D determined by the two methods is a measure of the error in the 
determination of D. In the similar three-dimensional system studied in ref. 0, D(T) was fit 
to the Vogel-Fulcher form 

D(T) = Ae- B ^ T - To \ (6) 

and we look for similar behavior. From the semilog plot of D versus 1/T shown in Fig. |3], 
it is easy to verify that the best fits are for T > 0. The best fit for all eight temperatures 
in Table |I| occurs for T ~ 1.8 and T ~ 1.7 using the values of D obtained from R 2 (t) and 
1/Q v (t), respectively. However, as shown in Table [TV] the results for T are sensitive to the 
number of data points which are included in the least squares fits. The range of reasonable 
fits illustrates the difficulty of obtaining meaningful values of To. We conclude that the 
temperature-dependence of D(T) is consistent with the Vogel-Fulcher form with To in the 
range 1 < T < 2. Note that the Vogel-Fulcher form of -D(T) implies that the system loses 
ergodicity at T = T . 

We also calculate the single particle energy fluctuation metric Q e (t) (see ([!])) to see if it 
exhibits ergodic behavior. If the system is ergodic, we expect that |10[ 



Q e {t) -> 2r e /t. (t » 1) (7) 

We interpret r e as the energy mixing time. The energy fluctuation metric was computed 
"on the fly" in the same way as fl v (t), but for a duration of 50 000 r rather than 100 000 r. 
We observed that l/Q € (t) becomes approximately linear even at the lowest temperature 
T = 2.25, indicating that the system is ergodic according to this measure. Our results 



for the T-dependence of t £ (T) based on (|7]) are summarized in Table |l| and are plotted 



in Fig. f|. Note that the largest mixing time is order 10 000 r. For comparison, we also 



show in Table [III] the values of td = o j (4D), the mean time it takes a particle to traverse 
the distance a. As summarized in Table [TV], the fits of r e (T) to the Vogel-Fulcher form 
T e (T) = Ae~ B ^ T ~ T ^ in different temperature intervals also yields a range of values for the 
parameter T e . We conclude that the temperature-dependence of r e is consistent with the 
Vogel-Fulcher form with a value of T t consistent with the value of T determined from the 
self- diffusion coefficient. 

3. Evidence of pseudospinodal behavior 

Now that we have established that the two component, two-dimensional supercooled 
Lennard- Jones system has properties similar to those observed in other simulations of deeply 
supercooled systems, we explore how the mean-field interpretation discussed in Section 1 is 
applicable to the present Lennard- Jones system. If the Lennard- Jones system exhibits pseu- 
dospinodal effects, we should find behavior analogous to that observed in the mean-field glass 
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model || and in Ising models with long, but finite range interactions 0. In these systems, 
the static structure function S(k) appears to diverge at a nonzero value of k if its behavior 
is extrapolated from high T or small magnetic field, but the extrapolated singularity is not 
observed if measurements are made too close to the apparent singularity. 

We compute S(k) (for the majority particles) from the saved particle positions using the 
relation 

1 N 

SW = ]d£e ik - r f- (8) 

i=l 

Angular averages were computed by considering all possible fc-vectors. We calculated S(k) 
for the same configurations as were used for obtaining g(r). The fc-dependence of S(k) 
at T = 2.7 is plotted in Fig. ^. Because the number of particles is fixed, we expect that 
S(k) — > 1 as k — > 0. The important feature of S(k) is the diffraction peak at k = k « 7.15. 
We interpret the height of this peak as a fc-dependent susceptibility x(fc ,T), and the width 
w(ko,T) as an inverse length proportional to the size of the correlated regions in the liquid. 
These quantities were extracted from S(k) by fitting the region around k = ko to the four 
parameter form, S(k) — a + b/[(k — k ) 2 + c]. 

Our results for x{ko,T) and w(ko,T) are listed in Table | and plotted in Fig. |6|. Note 
that x(ko,T) increases by a factor of approximately 1.6 and w(ko,T) decreases by a factor 
of approximately 2.0 as T is lowered from 5.5 to 2.5. Because this range of T is limited, we 
can fit x{ko,T) and w(ko,T) to a variety of functional forms. Given the divergent behavior 
of the first peak of S(k) in the mean-field model discussed in Section 1, we look for fits 
to the functional forms, x(fc ,T) ~ (T — T s )~ 7 and w(k ,T) ~ (T — T s ) u , where T s , 7, 
and v are fit parameters. We find that the best fit to the assumed power law form is 
X(k , T) ~ (T - 2.6)~ ai6 and w(k , T) ~ (T - 2.6) a29 if results in the interval 2.7 < T < 4.6 
are included. A least squares fit in the interval 3.1 < 5.5 yields x(ko,T) ~ (T — 3.0) -0 08 and 
w(ko,T) ~ (T — 3.0) 0,19 . Given the limited range of T, these power laws fits are justified 
only in the context of our rigorous results for S(k) in the mean-field model for which 
x(k ,T) ~ (T — Ts) -1 . The power law fits suggest that the increase of the height and the 
decrease in the width of the first peak of S(k) is influenced by a weak singularity at T = T s 
with T s in the range 2.6 <T S < 3.0. As expected, no evidence of a singularity is found if we 
fit the results for x{ko,T) in the interval 2.25 < T < 5.5. 

4. Cluster scaling and lifetime. 

Given the preliminary pseudospinodal interpretation that we presented in Section 3, we seek 
other evidence for the existence of a pseudospinodal and mean-field behavior. For tempera- 
tures near the pseudospinodal the system should show signs of an instability and the system 
should partially phase separate. That is, we should see regions where the majority particles 
dominate and locally order. For this reason we look for long-lived structures whose con- 
stituent particles remain in close proximity to each other over extended times at sufficiently 
low T. Because of the strong repulsive Lennard- Jones interaction near the origin, these 
structures will not be identical to the clumps found in the mean-field model, but instead 
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should be caused by a cage effect. A visual examination of the configurations shows evi- 
dence of a partial phase separation in which a significant fraction of the majority particles 
form clusters of triangular-like structures which become better defined as the temperature 
is decreased. Qualitatively, the lifetime of these visual clusters in the range of temperatures 
of interest appears to be much longer than the period of oscillation of the particles about 
their mean position in the clusters. 

To find these structures we seek to define a "solid-like" particle such that clusters of these 
particles have the above qualitative properties. However, unlike the Ising model near the 
critical point [H| and near the spinodal [T5|], there is no theoretical definition of clusters in 



a dense liquid, and we have to rely on our physical intuition to define them. We will assume 
that a cluster consists of a group of solid-like particles and that if two solid-like particles 
are near neighbors, they belong to the same cluster. Ideally, we would like to introduce a 
physical property that exhibits bimodal behavior with one peak corresponding to solid-like 
particles. However, all of the physical quantities we measured exhibit a single peak, and 
hence we will need to introduce a cutoff to distinguish solid-like and non-solid-like particles. 
Nevertheless, we will find that the properties of the resultant clusters do not depend strongly 
on the choice of cutoff. 

Because we are interested in the local environment of the particles, it is natural to do a 
Voronoi construction [T7J and determine the Voronoi polygon of each particle and its Voronoi 



neighbors. Quantities associated with the Voronoi construction include the distribution of 
the number of edges of the Voronoi polygons, the distribution of the area of the polygons, 
and the distribution of the length of the sides. For our two-dimensional system, the mean 
number of edges of the Voronoi polygons is six and the distribution of the number of edges is 
temperature-independent. This independence is expected because the Voronoi construction 
in two dimensions is insensitive to thermal fluctuations |TB| . 



To quantify the temperature-dependent changes in the distribution of lengths of the sides 
of the Voronoi polygons with six sides, we introduce the standard deviation at of the edge 



length of a particular particle as [pl| 



of = (?) - (V, (9) 

where (f(i)) = (1/6) J2a=i fi^a) and £ a is the length of edge a of the hexagon of interest. 
If a particle were in a perfect triangular environment, then at = 0. In Fig. |7| we show the 
estimated probability density P{h) of the quantity h = a^/(£) 2 , a measure of the "hexag- 
onality" of the polygons. At T = 5.5, the Voronoi polygons are irregular, giving rise to a 
high percentage of particles with h > 0.1. At T = 2.7, the most probable value of h is at 
h = 0.014, and many more particles have regular Voronoi polygons. 

In the following, we define a majority particle to be solid-like if it has six Voronoi neigh- 
bors and if the condition h < 0.1 is satisfied. The choice of the cutoff parameter is not crucial 
and the qualitative properties of the clusters are independent of the cutoffs over a wide range 
of values |18|| . A typical configuration of the system at T = 3.3 is shown in Fig. || 

Given our criteria for solid-like particles and our definition of clusters, we can determine 
the properties of the clusters. A separate program using a Voronoi construction and a 
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Hoshen-Kopelman cluster labeling algorithm was used to determine the clusters. Because 
the height of the maximum of S(k) increases and the width of the peak decreases as T is 
lowered, we expect that the mean size of the clusters grows as T is lowered until their growth 
becomes "frustrated" by the presence of the larger minority component and by the different 
orientations of the other clusters. (The size of a cluster is its number of particles.) The 
behavior of the cluster size distribution n s as a function of the size s is plotted in Fig. |9| for 
T = 3.3 and N = 500. The cluster distribution is normalized by the number of (majority) 
particles so that n s is the probability that a particle belongs to a cluster of size s. The results 
for n s are averaged over 100 000 r. We expect that the presence of a pseudospinodal leads to 
power law scaling of the clusters if T is close, but not too close to the pseudospinodal, that 
is, we expect n s ~ s~ x for a range of values of s. The lack of a true spinodal should lead to 
a cutoff or lack of scaling for large clusters or for temperatures far from the pseudospinodal. 
From the log-log plot of n s versus s in Fig. |9| we see that the s-dependence of n s is consistent 
with a power law dependence with an exponent of x ~ 1.8 over approximately one decade 
of s. 

The behavior of n s for lower temperatures is qualitatively different. For example, compare 
the behavior of n s averaged over the first 50 000 r of the run at T = 3.1 to the plot of n s 
averaged over the second 50 000 r of the run (see Fig. |TD|) . The difference in the cluster 
size distribution for larger values of s indicates that the lifetime of the clusters at T = 3.1 
is the same order of magnitude as the duration of our runs. The fact that larger clusters 
require a longer time to reach equilibrium than smaller clusters has been found previously in 
temperature quenches of Ising models near the spinodal |T3| . Our simulation results for n s 
for lower values of T show similar behavior. We will discuss other estimates of the lifetime 
of the clusters in the following. 

The behavior of n s for T = 5.5 is shown in Fig. |TT|. No evidence for simple power law 
behavior is observed, but the s-dependence of n s is consistent with a fit to the assumed 
form, n s oc s ~ x e~ s ^ ma with x — 1.1, and m s = 13. Similar fits can be made for T = 4.7 with 
x = 1.2, and m s = 21, and T = 3.7 with x = 1.4, and m s = 42. We note that the effective 
value of the power law exponent x and the cutoff parameter m s increases as the temperature 
is decreased, suggesting that a weak singularity is being approached. 

Because our results for n s for N = 500 might be affected by finite size effects, we did 
a run for N = 20 000 particles at T = 2.7 for a time of 30 000 r. The system was run for 
15 000 t before data was collected. A log-log plot of n s versus s is shown in Fig. 0. It is 
clear that the values of n s for larger s are not in equilibrium. However, a power-law fit of 
n s in the range 4 < s < 61 yields a slope of x — 1.75, a value of x consistent with the value 
estimated from our results for N = 500 at T = 3.3. 

We expect that if the clusters are important near the glass transition, their lifetime should 
increase as the glass transition is approached. As indicated in Fig. [H| we know qualitatively 
that the lifetime of the larger clusters becomes very long as the temperature is reduced. We 
introduce a measure of the cluster lifetime by measuring the time-dependence of a metric 
associated with the solid-like particles. We divide the simulation cell into boxes and compute 
the number of solid-like particles in each box. The idea is to find if the time averaged number 
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of solid-like particles in each box becomes the same when the time t » 1. We take n a to be 
the number of solid-like particles in box a and compute the cluster fluctuation metric fl c \(t) 
defined as 



1 



^(t) = — J2[n a (t)-(n(t))Y- (10) 



In (|Tol) n a (t) is the mean number of solid-like particles in box a at time t, and (n(t)} is the 
mean occupancy averaged over all Nf, boxes. For our runs we divided the system into 5x5 
boxes. A linear increase in 1/Q c \(t) defines the cluster mixing time r c \ as in (|7|). 

At T = 5.5, Q c \(t) exhibits ergodic behavior with t c \ ~ 500 r (see Fig. |13|a). Our estimates 



for r c i for our runs are summarized in Table |T|. Note that for T = 3.1 (see Fig. |T3|b), r c i is 
order 5 x 10 4 r, a time which is comparable to the duration of runs. For T < 3.1, our estimates 
of r c i are longer than the duration of our runs and are not meaningful. We interpret the time 
r c i as an estimate of the lifetime of the clusters. Note that at T = 3.1, the times associated 
with the cluster lifetime and the motion of single particles differ by a factor of 10 3 . The In- 
dependence of r c i is best approximated by a Vogel-Fulcher form r c \ ~ e c ^ T ~ Tcl \ where T c \ is 
the extrapolated temperature at which the cluster lifetime would become infinite. Although 
our estimates for r cl are only qualitative, the estimated value of T cl found by considering the 



values of t& in the range 3.1 <T< 5.5 yields a reasonable fit with T& ~ 1.6 (see Fig. [14]). 



This estimate of T c \ does not vary much if the results at T = 3.3 and T = 3.7 are omitted 
(see Table |TV|). 

Another single particle decorrelation time can be extracted from nf,(t), the number of 
unbroken Voronoi bonds remaining after a time t. At t = 0, only Voronoi bonds between 
small particles that have exactly six small neighbors are counted. If at a time t later, there is 
no longer a bond which joins the same pair of particles, the number of bonds is reduced, and 
we do not consider this pair of particles again. We find that (rib(t)nb(0)) ~ e~ l l Th . Because 
a bond is broken every time a particle changes neighbors, we expect that r& ~ td- We 
computed Tf, for only a few temperatures and found that Tf, is comparable to tjj. 

Given the very long lifetime of the clusters for T < 3.1, it is difficult to make estimates of 
the errors associated with various quantities. For example, even though we made long runs 
at low temperatures, we found only one statistically independent configuration as far as the 
clusters are concerned. For this reason, quantities which are measures of the structure of 
the system, such as S(k), are probably not adequately sampled for T < 3.1. And although 
our measures of single particle properties such as the velocity metric show the system to be 
ergodic at this level, the values of D at low T might also be inadequately sampled because 
the motion of the particles is influenced by the presence of the long-lived clusters. 

It is not clear from our results and various fits whether we can identify one or two temper- 
atures which can be associated with a glass transition. Our estimates for the temperature T 
at which the self-diffusion coefficient vanishes and measures of the single particle properties 
become nonergodic are in the range 1 < T < 2. In comparison, our estimates for the tem- 
perature T s at which the fc-dependent susceptibility x(^o) an d the inverse width w~ 1 (/c ,T) 
would diverge if a spinodal were present are in the range 2.6 < T s < 3.0. The cluster distri- 
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bution exhibits power law scaling at a temperature which is close, but not too close, to our 
estimate of T s . At T pa 3.1, the lifetime of our clusters is already comparable to the duration 
of our runs. Nevertheless, if we extrapolate the cluster lifetime to lower temperatures by 
fitting the cluster lifetime to a Vogel-Fulcher form, we find that the cluster lifetime becomes 
infinite at T c \ ~ 1.6, a temperature consistent with our estimate of To. 

Given the small size of our system for all but one of our runs and the limited number 
of temperatures available, it is difficult to make quantitative conclusions in spite of the 
relatively long duration of our runs. Moreover, as we have emphasized, our interpretation of 
our molecular dynamics data is justified only in the context of our rigorous results for the 
mean-field model of a structural glass discussed in Section 1. 

5. Mean-field interpretation and summary. 

On the basis of our molecular dynamics simulations, we can conclude that the deeply su- 
percooled, two component, two-dimensional Lennard- Jones system is heterogeneous. The 
heterogeneity can be characterized both dynamically and statically, that is, there are long- 
lived spatially correlated regions of all sizes up to a cutoff. This qualitative conclusion is 
consistent with the results of rotational diffusion experiments on probe molecules in su- 
percooled o-terphenyl [fJTJ, which indicate that the dynamics in supercooled o-terphenyl is 



spatially heterogeneous, and with the results of other simulations [ POU . 

An important question is the appropriate theoretical understanding of the origin of the 
spatial heterogeneity. In a future publication we will give a mean-field argument for the 
origin of the scaling behavior of the cluster size distribution in terms of an arrested nucleation 



picture 22 



We have found indirect evidence for mean-field behavior and the influence of a pseu- 
dospinodal. As discussed in Section 2, the height of the diffraction peak of the static struc- 
ture function S(k) for the mean-field structural glass model exhibits a true divergence Hj in 
the limit R — > oo. For finite range R, the peak of S(k) appears to diverge if it is measured for 
values of the temperature T which are close, but not too close, to the apparent singularity 
and the data is extrapolated to lower T. However, a singularity in the peak of S(k) is not 
observed if measurements are made too close to the apparent singularity. This behavior of 
S(k) is characteristic of a pseudospinodal. Is there a spinodal in the Lennard- Jones sys- 
tem? The answer is no, because the range of the Lennard- Jones potential is finite. However, 
we found that the height and the inverse width of the diffraction peak of S(k) exhibits a 
weak power law divergence if their behavior is extrapolated from high T. This behavior is 
consistent with a pseudospinodal interpretation. 

We expect that the Lennard- Jones system would be better described by mean-field the- 
ory as the density is increased, because the number of interactions each particle experiences 
increases. At present, we do not know how to calculate the effects of the pseudospinodal in 
dense Lennard- Jones systems ||23|| , and we need to rely on simulations. Several other mea- 



surements suggest that the deeply supercooled, dense Lennard- Jones liquid can be described 
at least qualitatively by a mean-field picture. Glaser and Clark (T^| found cluster scaling in 



a simulation of a two-dimensional Lennard- Jones system near the freezing transition. On 
the basis of a mean-field theory, it has been predicted that near the liquid-solid spinodal, 
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the nucleating droplets are fractal-like rather than compact objects [p4p5| . This effect has 



been observed in simulations of a three-dimensional Lennard- Jones system [ 26fl . A third 



result consistent with a mean-field interpretation is that the measured value of the fractal 
dimension of the structures formed in a single component two-dimensional Lennard- Jones 
system undergoing spinodal decomposition is consistent with that predicted by mean-field 
theory [ 15|,27| . These results, together with the results reported here, suggest that a mean- 
field interpretation is applicable to dense Lennard- Jones systems under the proper conditions. 

We note that although the clumps in the mean-field glass model and the solid-like clus- 
ters in the Lennard- Jones system have some properties in common, for example, their long 
lifetime, the clumps are not directly analogous to the clusters. As discussed in ref. the 
clump size distribution is a Gaussian in contrast to the power law distribution that we found 
for the solid-like clusters. Although a theoretical definition of the clusters in the mean-field 
glass model does not exist, we can follow a similar approach and assume that a cluster in the 
latter system is a group of clumps such that each clump is in a triangular environment (in 
two dimensions). A preliminary investigation |28| of such a criterion for a cluster of clumps 
yields clusters which have a power-law distribution near the spinodal consistent with our 
results for the Lennard- Jones system. 

Our argument for the association of the glass transition with an underlying thermody- 
namic transition is consistent with the recent interpretation by Nagel and coworkers |^9| 



of the frequency dependence of e"(u), the imaginary part of the dielectric susceptibility, in 
organic glass forming liquids. Nagel and coworkers [29] have fitted e"{uj) to a single scaling 
curve over 13 decades of frequency for a wide range of T and for many glass formers. If the 
temperature-dependence of the high frequency, power law behavior of e"(uj) is extrapolated 
to lower frequencies, they found that the static susceptibility diverges at a temperature T a 
with the corresponding mean- field exponent JJtJ. The magnetic susceptibility in a dipolar- 
coupled Ising spin glasses is found to have similar behavior. Given that these experiments 
have been done on fragile glass forming liquids with long molecules and on dipolar magnets, 
a mean-field interpretation of these results is consistent with our picture of the glass tran- 
sition. That is, these systems are well approximated by mean-field models due to the large 
number of simultaneous interactions which each molecule has with its neighbors f31 |. 



Our interpretation of our simulations in terms of an underlying thermodynamic transition 
in Lennard- Jones systems is consistent with the results of recent laboratory experiments. 
However, the interpretation of the single particle behavior in terms of a distinct kinetic 
transition is more open to question. Although such an interpretation is rigorous for our mean- 
field model of a structural glass, Nagel and coworkers |30| have interpreted their experimental 
results in terms of a single glass transition temperature. That is, the temperature T a at which 
the static susceptibility is extrapolated to diverge is the same temperature at which the self- 
diffusion coefficient is extrapolated to vanish. In contrast, we find two distinct temperatures 
in the mean-field glass model and also can interpret our simulation results of the Lennard- 
Jones system in terms of two temperatures. We note that in the mean-field limit the number 
of particles in each clump is infinite, and the system would remain indefinitely in a local free 
energy minimum as determined by the number and location of the clumps. The duration 
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of our simulations of the Lennard- Jones system is sufficiently short (~ 10~ 7 s) that at low 
temperatures the clusters do not diffuse and particles in the interior of the clusters are 
trapped. Moreover, the mean cluster lifetime near the glass transition is estimated to be an 
order of magnitude longer than our longest runs. This picture of static clusters is consistent 
with the mean-field model, but might not be appropriate for experimental time scales. Hence, 
it is possible that if we were able to run for longer times, we would find that the extrapolated 
temperature at which the self-diffusion coefficient vanishes and the extrapolated temperature 
of the underlying thermodynamic glass transition would approach each other. 

Other workers have also interpreted the behavior of supercooled liquids near the glass 
transition in terms of clusters. Kivelson and coworkers |32] have proposed a thermody- 
namic theory of supercooled liquids based on the assumed existence of a "narrowly avoided" 
thermodynamic phase transition. The avoided transition is attributed to the existence of 
strain. In contrast to this transition, the spinodal transition always occurs below the first- 
order transition temperature. The reason that the spinodal transition is "avoided" in our 
interpretation is due to the fact that the system is not really mean-field. 

On the basis of our results for the temperature dependence of w(k,T), the width of 
the first peak of S(k), we can conclude that there is an increasing length scale which is 
associated with the clusters as the temperature is decreased. In our interpretation this 
increasing length scale is due to the effects of the pseudospinodal. This length scale would 
diverge if a spinodal were really present. The relation of this increasing length scale to the 
increasing maximum length scale for propagating transverse current correlations observed 
by Mountain |33j requires further study. In contrast, an increasing correlation length was 



not observed in a simulation [51] of the translational and orientational correlation functions 
of a two-component three dimensional Lennard- Jones system. 

The effects of the pseudospinodal and the incipient thermodynamic glass transition will 
be more or less apparent depending on the interaction range, the details of the interaction, 
and the spatial dimension |7[ . We do not expect to find spinodal-like effects in all supercooled 
liquids. These considerations suggest that there is a class of materials for which the observed 
glass transition is associated with a pseudospinodal and an incipient thermodynamic tran- 
sition, and other materials for which the observed glass transition might not be associated 
with such effects. 

Based on our Monte Carlo and theoretical studies of a mean-field model and our molecular 
dynamics results for a two-dimensional Lennard- Jones system, we suggest that the latter is 
in the class of systems whose behavior can be attributed to an incipient thermodynamic 
instability (the pseudospinodal). We emphasize that a true thermodynamic glass transition 
does not exist in the Lennard- Jones system, even though the pseudospinodal has measurable 
effects including increasing length and time scales as the pseudospinodal is approached. In 
addition to this glass/pseudospinodal transition, there is a temperature (for fixed density) 
which can be interpreted as a kinetic transition below which the self-diffusion coefficient is 
not measurable during our observation time. 

We are presently simulating much larger Lennard- Jones systems in two dimensions to 
obtain better statistics for the clusters over a wider range of sizes and over a range of 
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temperatures above the glass transition. If our mean-field interpretation is correct, we should 
be able to observe similar behavior in three dimensions where mean-field behavior should be 
even more apparent. However, the identification of the clusters in three dimensions is not 
straightforward because of the existence of a variety of possible local symmetries which are 
more affected by thermal fluctuations than in two dimensions. 
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Appendix A. The enthalpy. 

Many computer simulations of glasses show that the heat capacity Cp has a maximum in the 
vicinity of the glass transition. For example, Wahnstrom has computed the temperature- 
dependence of the energy of a two-component, three dimensional Lennard- Jones system and 
has found that the slope changes at a temperature where the dynamical anomalies are most 
pronounced. In the following, we discuss our measured values of the temperature dependence 
of the enthalpy H. 

We measured H using constant pressure molecular dynamics [|3£] for iV = 500 particles. 



The system was equilibrated at T « 10 and a series of measurements of H were performed 
at progressively lower temperatures spaced approximately 0.02 apart for the higher values of 
T and approximately 0.01 apart at the lower values of T. At each value of T, the system was 
equilibrated for 250 r and H was averaged over the following 250 r. These runs are relatively 
short in comparison to the runs reported in the main text. 

Our results at P = 70 for H and Cp are shown in Fig. [15]. A careful inspection of H(T) 
shows that its slope changes as a function of T. Note that Cp(T) increases as T is lowered, 
reaches a maximum at T w 4 and decreases as T is lowered further. The slope of H(T) was 
computed as each value of T = Tj from the numerical derivative, [H{Ti + i)—H{Ti)]/\Ti + i—Tj\. 



The values of CpiT = Tj) shown in Fig. [15| were computed by doing a least squares fit to 15 
successive slopes in the interval T i+7 — Tj_ 7 . If we consider less than 15 points, the results 
for Cp(T) were too noisy. The results for more than 15 points tended to smooth the peak 
in Cp. 

We note that the spinodal singularity in the mean-field glass model is at k = ko ^ 0, 
and hence we do not expect thermodynamic quantities such as Cp to exhibit a divergence 
at the spinodal in the mean- field limit. However, for a system that does not exhibit a true 
spinodal such as the present Lennard- Jones system, it is possible that the coupling between 
the k = k and k = modes might lead to a maximum in quantities such as Cp near the 
pseudospinodal. 
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TABLES 



TABLE I. Values of the mean temperature T, the mean pressure P, the height of the first peak 
of the radial distribution function g(r), and the height x(^o) an d width uu(ko) of the diffraction 
peak of the static structure function S(k) averaged over runs of 100 000 r each. The mean pressure 
was computed from the virial. The width of the first peak of g(r) does not appear to change with T 
and is not listed. The density p is fixed at p = 0.952. As stated in the text, the lowest temperature 
run is probably not in complete thermal equilibrium. 



T 


P 


S^max) 


X(fco) 


w(k ) 


5.5 


91.1 


4.30 


2.59 


1.11 


4.6 


85.3 


4.54 


2.70 


0.95 


3.7 


77.8 


4.91 


2.87 


0.87 


3.3 


74.9 


5.15 


3.10 


0.76 


3.1 


72.9 


5.33 


3.49 


0.63 


2.7 


69.0 


5.68 


4.28 


0.45 


2.5 


67.3 


5.86 


3.94 


0.55 


2.25 


65.0 


6.07 


3.54 


0.54 



TABLE II. Summary of results for the self-diffusion coefficient D as a function of temperature 
T at constant density p = 0.952. The second column represents the estimates of D determined 
from the slope of R 2 (t), and the third column gives the estimates of D from the slope of 1/Q v (t) 
(see (§)). 



T 


D [from R 2 (t)] 


D [from l/n v (t)} 


5.5 


4.4 x 10~ 2 


4.8 x 10- 2 


4.6 


3.0 x 10~ 2 


3.1 x 10~ 2 


3.7 


1.2 x 10~ 2 


1.2 x 10~ 2 


3.3 


6.9 x 10~ 3 


6.8 x 10~ 3 


3.1 


4.1 x 10~ 3 


4.2 x 10~ 3 


2.7 


1.1 x 10~ 3 


8.1 x 10~ 4 


2.5 


4.5 x 10~ 4 


4.2 x 10~ 4 


2.25 


1.1 x 10~ 4 


6.7 x 10~ 5 
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TABLE III. Summary of characteristic times. The cluster mixing time r c i is computed as in 
(|7|) from the cluster metric f2 c i. The single particle time td = cr2 /(4-D), the mean time it takes a 
particle to diffuse a distance a, and the energy mixing time t £ from (0) are shown for comparison. 
(The values of td are obtained from the second column in Table |TJ.) 



T 


Tel 


TD 




5.5 


5 x 10 2 


5.7 


16 


4.6 


1 x 10 3 


8.4 


28 


3.7 


5 x 10 3 


20 


69 


3.3 


2 x 10 4 


36 


140 


3.1 


5 x 10 4 


60 


260 


2.7 


5 x 10 5 


230 


1100 


2.5 


5 x 10 5 


560 


1400 


2.25 




2300 


11000 



TABLE IV. Range of fits for the temperature parameter To from least squares fits of the 
self-diffusion coefficient D, the energy mixing time r e , and the cluster mixing time T c \ to the 
Vogel-Fulcher form (|6|). 



temperature interval 


T [from R 2 (t)] 


T [from l/n v (t)} 


2.25 < T < 5.5 


1.82 


1.67 


2.7 < T < 5.5 


1.83 


1.21 


2.25 < T < 4.6 


0.97 


1.30 


2.7 < T < 4.6 


0.94 


1.30 


T e 


2.25 < T < 5.5 


no fit 




2.7 < T < 5.5 


1.25 




2.7 < T < 4.6 


1.25 




3.1 < T < 5.5 


1.68 




3.1 < T < 4.6 


1.95 




T d 


2.25 < T < 5.5 


no fit 




2.7 < T < 5.5 


1.58 




2.7 < T < 4.6 


1.57 




3.1 < T < 5.5 


1.58 
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FIGURES 



FIG. 1. Plot of the radial distribution function g(r) versus r at T = 2.7. Note the split second 
peak, a possible signature of a glassy state. Only the small, majority particles were included. The 
dimensionless quantity r is measured in terms of the Lennard- Jones parameter a. 



FIG. 2. Plot of 0, v (Q)/Q v (t), the reciprocal of the velocity fluctuation metric, at T = 2.7. Note 
that l/£l v (t) is not linear for t < 250 r, but becomes approximately linear for longer times. The 
temperature T is dimensionless (see text). 



FIG. 3. Semilog plot of the (dimensionless) self-diffusion coefficient D versus 1/T. The data 
points (filled circles) are taken from the second column in Table ||. The solid line represents the 
best fit to the results for D in the interval 2.25 <T< 5.5 extracted from R 2 (t) to the Vogel-Fulcher 
form D(T) = Ae~ B ^ T - T ^ with A = 0.16, B = 4.7, and T = 1.82. 



FIG. 4. Semilog plot of the (dimensionless) energy mixing time r e defined in (^) versus 1/T. 
The data points (filled circles) are taken from Table III. The solid line represents the best fit to the 
Vogel-Fulcher form r e = Ae B ^ T ~ T ^ with A = 1.24, B = 9.8, and T t = 1.25 for the six temperatures 
in the interval 2.7 < T < 5.5. 



FIG. 5. Plot of the static structure function S(k) at T = 2.7 as a function of k. The height of 
the first diffraction maximum increases and its width decreases as T is lowered from T = 5.5 (see 
Table IB). 



FIG. 6. (a) The temperature dependence of x(^o> T), the height of the diffraction peak of S(k), 
at k = ko ~ 7.15. The solid line represents the best fit in the range 2.7 < T < 5.5 and has the form 
(T — 2.6)~ ' 16 . (b) The temperature dependence of w(ko,T), the width of the diffraction peak of 
S(k). The solid line represents the best fit in the range 2.7 < T < 5.5 and has the form (T — 2.6) ' 29 . 



FIG. 7. Plot of P(h), the estimated probability density of h, the relative variance of the edge 
length of the Voronoi polygons, at T = 5.5 (solid line) and T = 2.7 (dotted line). The data points 
are not shown to avoid confusion. Regular hexagons and hence lower values of h are much more 
likely at low temperatures. Note that P(h) is not bimodal at any temperature. The shoulder in 
P(h) at h Ri 0.25 appears to be real rather than statistical error. 



FIG. 8. A typical configuration at T = 3.3 showing the Voronoi polygon for each particle and 
the clusters of the majority small particles (shaded hexagons). 
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FIG. 9. Log-log plot of n s , the cluster size distribution, versus size s at T = 3.3 for N = 500 
particles. The solid line with slope x = 1.85 is the best fit in the range 6 < s < 61. 

FIG. 10. Log-log plot of n s versus s at T = 3.1 for TV = 500 particles. The results for the first 
50 000 t (open circles) and the second 50 000 r (open squares) are shown separately. The system 
was run for 50 000 r before data was taken. Note the difference in the size distribution of clusters 
for larger values of s. 

FIG. 11. Log-log plot of n s versus s at T = 5.5 for N = 500 particles. The plot shows 
curvature indicating that there is not a simple scaling regime. The solid line is a fit to the form 
n s = As- x e~ s / ms with A = 0.038, x = 1.1, and m s = 13. 

FIG. 12. Log-log plot of n s versus s at T = 2.7 for N = 20 000 particles. The linear fit was 
done for 5 < s < 61 and yields a slope of x = 1.75. 

FIG. 13. The time-dependence of f2 c i(0)/fJ c i(i), where f2 c i(t) is the cluster fluctuation metric 
defined in (|T(i|) at (a) T = 5.5 and (b) T = 3.1. Note the very different vertical scales in (a) and 
(b). The time t is given in terms of r defined in the text. 

FIG. 14. Plot of the cluster mixing time t& (solid circles) versus T, where t c \ is extracted from 
the linear time-dependence of l/f2 c i. The solid line represents the best fit to the Vogel-Fulcher form 
r cl (T) = Ae c ^ T - T ^ with A = 19.8, C = 11.8, and T cl = 1.58 in the interval 3.1 < T < 5.5. 

FIG. 15. The temperature-dependence of the (a) (dimensionless) enthalpy H and (b) the heat 
capacity Cp at P = 70. Note that Cp has a small peak at T ~ 4. As explained in the text, the 
slope of H(T) was computed at each value of T = Tj and the value of Cp at T = Tj was computed 
by doing a least squares fit to 15 successive slopes in the interval Tj + 7 — Tj_7. 
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